THREE-DIMENSIONAL STRUCTURE OF AN ALPHA ACCRETION DISK 
Wlodzimierz Kluzniak and David Kita 

Physics Department, University of Wisconsin, Madison, WI 53706, USA 

ABSTRACT 

An analytic solution is presented to the three-dimensional problem of steady 
axisymmetric fluid flow through an accretion disk. The solution has been obtained 
through a systematic expansion in the small parameter e = H/R (the ratio of disk 
thickness to its radial dimension) of the equations of viscous hydrodynamics. The 
equation of state was assumed to be polytropic. For all values a < 0.685 of the 
viscosity parameter, we find significant backflow in the midplane of the disk occuring 
at all radii larger than a certain value; however, in the inner regions of the disk the 
fluid always flows toward the accreting object. The region of backflow is separated 
from the region of inflow by a surface flaring outwards from a circular locus of 
stagnation points situated in the midplane of the disk. 

1. Introduction 

Accretion flows occur in a variety of astrophysical situations, often they take the 
form of a disk (e.g., Frank, King & Raine 1992). The first solutions for accretion 
disk flows were constructed numerically by Prendergast & Burbidge (1968), while the 
first analytic solutions were obtained by Shakura & Sunyaev (1973). Since then, it 
has been the custom in analytic, but frequently also in numerical, work to discuss 
essentially one-dimensional solutions, i.e. to obtain the radial structure of the disk by 
considering equations averaged over the thickness of the disk and only then to obtain 
an approximate "vertical" structure by separately considering equations describing 
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hydrostatic equilibrium (and possibly radiation transfer) in the direction perpendicular 
to the plane of the disk. This is also true in discussions of quasi-spherical flows, such as 
in the celebrated accretion-dominated flows (ADF) (e.g Narayan 1996 and references 
therein)^. 

As already shown by Urpin (1984) in a remarkable paper, consideration of the 
vertical gradients of the stress tensor leads to a solution in which the flow direction in 
the midplane of the disk is opposite to that in the subsurface layers. The flow cannot 
be properly described by its height-averaged value, a point dramatically evident in 
the numerical work of Kley & Lin 1992 who enforced spurious circulation flows in the 
meridional plane by adopting "height-averaged" boundary conditions at the edges 
of their computational domain (nevertheless, they found and correctly identified 
an outflow in the disk midplane, reminiscent of Urpin's solution). Several recent 
numerical calculations (e.g. Rozyczka, Bodenheimer & Bell 1994; Igumenshchev, 
Chen and Abramowicz 1995) also exhibit flows which can best be described in terms 
of a tendency for backflow to occur in the midplane of the disk. We believe this 
effect is not thermal in origin, and to investigate the dynamics of the phenomenon we 
solve analytically the three-dimensional equations of disk accretion using a polytropic 
equation of state for the fluid. 

Urpin (1984) included thermal effects but made the simplification of zero net 
angular momentum flow in the disk (equivalently, his self-similar solution is valid 
asymptotically for large radii). In our work we chose the opposite route — we neglect 

1 In an important contribution Narayan & Yi 1995 go beyond the one- dimensional 
solutions by numerically constructing axisymmetric ADF solutions which factorize the 
three-dimensional equations, i.e., solutions of the type f(r,9) = R(r)Q(9). However, 
the solutions we present in this paper are not factorizable. 
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thermal effects, but include the inner boundary condition — this allows us to exhibit 
the global character of the solution. In particular, we show how the backffow is fed by 
the inflowing fluid. In Section 2 we present the equations, in Section 3 we solve them. 
A discussion of the results is begun in Section 4 and concluded, in Section 5, with a 
detailed presentation of the velocity field in the disk. 
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1. Disk equations. 

1.1. Assumptions. 

We use cylindrical coordinates (r, 0, z) centered on the accreting object and make 
the following standard assumptions: 

(i) the gravitational force on a fluid element is characterized by the Newtonian 
potential of a point mass, 

*{r,z) = --j^, (1.1) 
with G the gravitational constant and M* the mass of the central star; 

(ii) the structure of the disk is symmetric under reflection about the (z = 0) 
midplane; 

(iii) the disk is in a steady state (d/dt = 0); 

(iv) the disk is axisymmetric (<9/<90 = 0), hence all quantities will be expressed in 
terms of the coordinates (r, z); 

( v ) M > kl; 

(vi) The disk is geometrically thin, i.e. \z\ r; 

(vii) Viscous torques are a small perturbation in the radial (r) and vertical (z) 
components of the equations of motion. 

Assumption i) implies also that the disk is not self-gravitating. The assumptions 
iii)— v) are consistent with the statement that the accretion time scale is much greater 
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than the Keplerian period. Assumption vi) implies that the rotational velocity is 
much greater than the local sound speed in the outer parts of the disk, ^> c s , and 
that the radial velocity is larger than the vertical one, |t> r | > \v z \. Assumptions v)-vii) 
taken together signify that the disk is approximately in hydrostatic equilibrium. 

Throughout this paper we will also assume that 

viii) the equation of state for the disk is that of a polytrope, i.e. 

P = Kp 1+1/n , (1.2) 

with n and K constant. Except in the Appendix we will take the polytropic index to 
be n = 3/2. 

For the inner boundary condition, we take vanishing of the viscous torque at 
some radius r m , corresponding to a maximum of the angular frequency, Q m , at 
the same radius. This boundary condition, which introduces into the problem a 
natural lengthscale, r + = Q^ n r^ l /(GM m ), is appropriate for black-hole disks and for 
stellar accretion disks about stars which are spinning-up (whether magnetized or 
not). For stars which are not spinning up, i.e. ones which transfer their angular 
momentum to the disk (such as non- magnetized stars rotating close to the equatorial 
mass-shedding limit and, possibly, for X-ray pulsars [accreting neutron stars] in 
their spin-down phase), the solution presented below is valid with the substitution 
1 — \Jr + /r — > 1 + yV + /r. Thus, our solution is universally valid for any thin, 
Keplerian accretion disk described by a polytropic equation of state. We expect that 
the qualitative features of our solution for the accretion flow will hold also for other 
equations of state. 

If the polytropic disk were in exact hydrostatic equilibrium, the angular frequency 
^ = v^/r would be constant on cylinders and it would be very easy to solve the 
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equations of motion, at least far from the inner boundary. In reality, viscous terms 
(which are of order o? in the alpha disk) break the hydrostatic equilibrium and cause 
the equations of motion to form a system of nonlinear, coupled, second order partial 
differential equations which are rather challenging to solve (even numerically), but 
which bring the reward of a solution whose salient features cannot be described by 
height-integrated, i.e. ordinary, differential equations. 



1.2. Equations of motion. 

We use the generalized Navier Stokes equations, along with the equation of 
continuity, to describe the accretion flow and to represent viscous interactions: 

P^+P^y- V)V^ = -VP-pV^ + V-a, (1.3) 

^ + V-(pV)=0, (1.4) 

where p is mass density, P is pressure, V is the velocity vector of a fluid element, and 
ip is the gravitational potential. The rank two viscous stress tensor, er, is assumed to 
have the following Cartesian components (Landau & Lifshitz 1959): 

'dVj dV k 2 



CTjk = V 



^6 jk V ■ V 



+ & jk V-V, (1.5) 



dxk dxj 3 

where £ is bulk viscosity and rj = up is the dynamic viscosity coefficient, both of which 
are functions of the coordinates. 

With the assumptions described in § |1 . 1| , the equations of motion in cylindrical 
coordinates become 

+ -1ft = -£-i^ + iF r , (1.6) 

or oz or p or p 
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where ip(r,z) is the gravitational potential given by eq. fll.lj ). F r and F 2 are 
respectively the r and 2 components of the divergence of the viscous stress tensor, i.e. 
the viscous force, and are given by: 
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i.n: 



1.3. Constants of integration. 



Vertical integration of eq. ( |1.9D , with the assumption of a steady state, yields 
an expression of the conservation of mass flow through cylinders. Usually this is 
written as M = — 27rrElv where M is the constant mass accretion rate through any 
cylinder (and hence onto the star), £ is the surface density in the disk, and tv is an 
effective (i.e. density-weighted, height-averaged) radial velocity. However, since we 
are interested in the z dependence of the radial velocity, v r , we choose to write this 
important equation as 



M 



-2nr 



pv r dz 



constant, 



1.12) 



where by convention M > for accretion, i.e. for v r < 0. The quantity M will serve 
as an integral of the motion for our accretion flow. 
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Another constant is obtained if, in the same spirit, we vertically integrate the 
angular momentum equation ( [L.7|) . If we first multiply both sides by r 3 and integrate 
over z from — oo to +00, we obtain: 



(rpV r )^^-dz- 



r^^ldz = l 
oz or 



_|_ r 3^ 



dz 



;i.i3) 



where in deriving the second term on the left hand side we have performed an 
integration by parts and set the boundary term to zero since p — > as z — > ±00. 

Using the equation of continuity (|1.9|) , we can transform the entire left-hand side 
into (d/dr) f_™ r 3 pv r fldz. The last boundary term involving dfl/dz also vanishes 
because rj = up ^ a,s \z\ — > 00 and, finally, integration over r gives 



J(r) - C 



-2-nr- 



V^-dz, 
or 



;i.i4) 



where —J(r) is the advection rate of angular momentum through a cylinder of radius 
r, C is a constant of integration, and the right-hand side is the net torque exerted 
by viscous interactions on the same cylinder. Note that this equation is exact for 
any azimuthally symmetric, steady flow in which no mass is exchanged through the 
surface at infinity (z = ±00), and in which no angular momentum is carried radially 
by radiation (Kluzniak 1987). 

Since we consider only cases when the accretion rate is never zero, we can 
introduce another constant j + = C/M. The torque vanishes when this new constant 
is equal to the height-averaged specific angular momentum j (weighted with radial 
momentum flux) or, correspondingly, when the height- averaged radial derivative of 
the angular momentum (weighted with dynamic viscosity) vanishes. That is, we can 
rewrite eq. (|1 . 14 ) as: 



M(j-j+) = -27rr 3 



rjdz 



dn 

dr ' 



;i.i5) 
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where j(r) = — 27rr 3 (f*™ pv r £ldz) = Mj, etc. If f2 is independent of z (i.e. constant 
on cylinders), the usual form of eq. ( |1.15|) is recovered by removing the bars. Thus 
j+ can be interpreted as the specific angular momentum at the zero-torque radius, 
r m . For this reason, we now define a new effective radius, r+, at which the Keplerian 
specific angular momentum is equal to j + , i.e. y/GM*r + = j + = f2(r m )r^. Note that 
in general the maximum value of Q is not equal to the corresponding Keplerian value, 
^( r m) 7^ Qkijm), and hence we do not expect r + to be the same as r m . 



1.4. The polytropic sound speed. 

In the standard theory of thin accretion disks, the local sound speed becomes of 
prime importance when modeling subsonic accretion. A clear advantage of employing 
a barytropic equation of state is that it reduces the number of variables by one. A 
polytropic equation of state also greatly simplifies calculation of the local sound speed, 
i.e. 

2 dP A l\P 



With the above relation we can rewrite the pressure gradients in eqs. (|1.6|) & (|1.8| ) in 
terms of c s , giving the following elegant expressions: 

ldP_J?c| \^L- n ^i ( 117 ) 

p dr dr p dz dz 

Now it is easy to show a basic result concerning a fluid in hydrostatic equilibrium 
in both the radial and vertical directions. Here eqs. ( |1.6| ) & ( |1.8| ), without the inertial 
and viscous terms, reduce to a simple form involving only c s , Q, and ip: 
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Taking d/dr of eq. ( 1.19|) and d/dz of eq. (|1.18|) , we obtain the familiar result that 
dQ/dz=0, i.e. Q is constant on cylinders for a polytrope in hydrostatic equilibrium 
(c/. Tassoul 1978). Since in eq. (|1.7|) the velocities are proportional to the viscosity, 
this already implies that in a barytropic disk of any thickness Q is independent of 
z to leading order in a Taylor expansion in the (small) viscosity parameter (except, 
possibly, when the specific angular momentum is constant, j = r 2 Q = const.). We will 
perform a systematic expansion in a different small parameter, the dimensionless disk 
thickness, but for subsonic flow the same zeroth order result will be recovered. 



1.5. Scaling the equations of motion. 

At this point it is of paramount importance that we scale all relevant quantities by 
their corresponding characteristic values. This will make the equations dimensionless 
and allow us to weigh the relative significance of each term that appears. Following 
Regev (1983) we scale all velocities (v z , v r , and c s ) by the characteristic sound speed, 
c s , all radial distances by some characteristic radius R (e.g. i?*), and all vertical 
distances by H, the typical vertical scale height in the disk. We also represent Q 
in units of (GM^/R 3 ) 1 ^ 2 = fifc*, the Keplerian angular velocity at the characteristic 
radius, and p in terms of a typical value p. Similarly, we scale the pressure by 
P — P Cs 2 , the kinematic viscosity by v = c s H, and the dynamic and bulk viscosity 
coefficient by ( = f] = v p. 

To apply a perturbative expansion technique to each equation we define an 
expansion parameter, e. Since we are interested in geometrically thin disks we choose 

e = £ = « 1, (1.20) 

R ^lk*R 

where we have used c s = HQk*, m agreement with the standard result from thin 
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disk theory that H ~ c s /Qk- in effect, e is a parameter which measures the relative 
"thinness" of the disk. 



Denoting the scaled forms of v r and v z by u and v respectively, we obtain the 
following set of non-dimensional equations: 

21 ^3/2 



2 du du ^ 2 

e u— — h ev— it r 

or dz 



e 2 (n 9 -^ 
\ dr 
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dr dz 
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r 
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pr dr 
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dv^ 
dr 



+ 



ld_ 

p dz 



p dz 



, 2 \ dv 



n d^ + l_d_ ( dv" 

dz pdz \ dz / 

e d ( du\ 



pr dr y dz 



e d . . d . . 



0. 



1.23) 



1.24) 



where we have used eq. (|1.17|) in rewriting the pressure gradients. Here eqs. ( |1.21|) , 
( |1.22|) , and ( |1.23|) are the scaled radial, angular, and vertical momentum equations 
respectively and eq. (|1.24j) is the scaled form of the continuity equation. Armed with 
the knowledge that e< 1 for a thin disk, we make eqs. ( |1.21|) - (|1.24|) the foundation of 
our analysis and proceed to perturbatively expand all dynamical quantities in powers 
of e. We will find u ~ e, v ~ e 2 , i.e., v r = 0(e 2 )v<p and v z = (9(e 3 )t>0; therefore the 
divergence terms V • V in eq. ( |1.5| ) contribute at order not lower than e 3 to eqs. 
and fll.8| ) — this formally justifies their frequent neglect. 
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2. Solution for the vertical structure by perturbative expansion in e. 



2.1. Power series in e = H / R. 



We expand all variables in powers of e and will evaluate eqs. (|1.21|) - (|1.24Q at 
various orders. We let 

(fi/fi*,) =fio + e^i + e 2 fi 2 + ..., (2.25) 

u = (v r /c s ) = u + eui + e 2 u 2 + (2.26) 

v = K/c~) = «o + e^i + e 2 ^2 + (2.27) 

as well as (c s /c s ) = c s0 + ec sl + e 2 c s2 + and (p/p) = po + epi + e 2 p 2 + with the 
assumption that the dimensionless vertical scale height of the disk, h(r) = H(r)/H, 
is of order unity, h ~ 0(1). All other variables like P, r), v, & M can be expressed in 
terms of these six fundamental quantities^. Our objective then is to calculate, order 
by order in e, the functional dependence of ft, u, v, c s , p, and h on the coordinates r 
and z alone. 



Assumption ii) of § regarding reflection symmetry about the (z = 0) midplane, 
implies that physical quantities such as Q, p, P, r], u, and c s are even functions of z, 
while v is odd under reflections through the equatorial plane. When we expand an 
even/odd function (e.g. Q) in powers of e <C 1, we require each term in the expansion 
(e.g. Qf, i = 0, 1,2,...) to be independently even/odd. 



2 In reality only five of these are independent for a polytrope since eq. ( |1.16 ) gives 
c s in terms of p (or vice- versa). 
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2.2. Viscosity-independent zeroth order results for the vertical structure. 

Examination of eq. ( |1.21 ) at zeroth order immediately gives 

Q = r~ 3 / 2 , (2.28) 



i.e. Qo is equal to the Keplerian value at the midplane, Qk = v^k/r = sjGM^/r 2, 
in conventional units. Though eq. ( |2.28j ) is consistent with the assumption of 
a rotationally supported disk, it cannot satisfy the inner boundary condition of 
n(r*) = f2* whenever the star rotates below its Keplerian value at the stellar radius, 
nor the more general zero-torque boundary condition. Clearly, our perturbative 
solution will be invalid in the limit r —* r + . This is because implicit in the scaling of 
§ p. . 5| has been the assumption that d/dr ~ e(d/dz). This approximation, however, is 



known to be patently false in the inner transition region between the central star and 
the Keplerian portion of the disk. 

Eq. ( |1.24| ), the equation of continuity, fixes v at zeroth order to be zero everywhere, 
v = 0. To see this more clearly, observe that d(p v )/dz = and so pov Q , the lowest 
order vertical component of the mass flow, is a function of r only. However, since v is 
odd with respect to the z coordinate, we know PqVq = for all points on the midplane 
(z = 0) and, not being a function of z, this product must then vanish everywhere. 
Clearly p ^ and thus v o = at all points in the rotationally supported disk. 

Moving on to first order in e for the angular velocity, we see that because Q is 
even with respect to reflections through the midplane, the first order correction to the 
angular velocity vanishes, Q± = 0. This result for f^i also has direct impact on the 
fluid velocities, since eqs. ( |1.22| ) & ( |1.24[ ) at order e now give u = v\ = 0. Using this 



result for u and v, we then find that the only surviving term of order e in eq. ( |1.23| ) 
involves the first order correction to the square of the sound speed, and thus (? sl = 0, 
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and hence p\ — P\ = 0, which is consistent with symmetry arguments for these three 
quantities. With this information, eq. ( |1.23|) simply becomes the standard equation of 
vertical hydrostatic equilibrium: 



1 dP d(c 



Po dz dz r 

We can solve eq. Q2.29Q , i.e. the vertical momentum equation up to corrections of 
order e 2 , and hence find c s0 , p , and P Q . In the case of polytropic index n = 3/2 we 
obtain the following relations (Hoshi, 1977): 



h 2 - z 2 

C *°M) = y —^3—, (2.30) 

(l2 2 \ 3/2 

. ( 2 - 31 ) 

//,2 _ 2\ 5/ 2 

fl»(r,z) = p8 /S =^^/-J • (2.32) 

Eqs. ( p.31| )- p.32j ) show that h(r) is now the height at which p = (and hence P = 0), 
implying that ft, is the true semi-thickness of the disk, though its functional dependence 
on r is still undetermined. The surface density, E(r), can also be derived in terms of h 
to lowest order in e: 

37T / ft 4 \ 

Eo(r) = ^, 0(iz = _^_j, (2 .33) 

where we have replaced the integration limits of ±oo by ±h, since by eq. ( |2.31| ) the 
polytropic disk terminates at z = ±h(r). 



It is now possible to put the integrated mass continuity equation (|1.12|) into 



dimensionless form. Scaling as before and using the fact that to lowest non-vanishing 



order in e, u = eu\ + ... and p = p + we can write eq. ( |1.12j) to lowest order in e: 

*+h r+h 



r+n r+n 

em = —r pudz ~ — er / p Q uidz, (2.34) 

J—h J —h 
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where fn = M/(2irpc s H 2 ) ~ 0(1), is a dimensionless constant. For an adiabatic 
index of n = 3/2, eq. (|1.2|) gives p ~ K~ 3 / 2 c s 3 and, since H ~ c s /Vtk*i we find 
m ~ M/c s 6 . Because m is independent of e for Mi 7^ 0, this implies that the unsealed 
mass accretion rate must scale as M ~ £*(e 6 ), and hence M depends sensitively on the 
relative "thinness" of the disk. As we shall soon see, the scaled eq. (|2.34|) is of prime 
importance in determining the vertical structure and must therefore be included with 
eqs. (imP-QOl. 



Up to this point, all the results obtained in this section, including eqs. ( |2.28|) and 
( |2.30D -( p.33; ), are common to the outer regions (r >> r + ) of any standard thin disk 
with polytropic equation of state, for any viscosity prescription. However, all of these 
expressions depend intrinsically on the vertical scale height, h(r), which cannot be 
determined as an explicit function of r without a form for i](r, z) being first specified. 
We cannot obtain solutions for the lowest non- vanishing orders in e of v r and v z , nor 
can we evaluate the first nonzero correction to Q without specifying the viscosity 
prescription. 



2.3. Lowest-order results using the standard a-disk prescription. 



Let us continue all calculations under the assumption that the viscosity is given 
by the a-disk formulation. In view of the large uncertainty in modeling the kinematic 
viscosity, authors in the past have generally neglected any z dependence for v. 
Laboratory studies of turbulent jets undergoing free expansion lend some support to 
this hypothesis (cf. Urpin 1984a, Monin & Yaglom 1965). However, we find this 
to be unacceptable when solving for v r (r, z) and v z (r,z) in a polytropic disk as it 
leads to divergent expressions at the surface of the disk, v (r, ±h) — > 00 (Kita 1995). 
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We choose then to modify the alpha prescription by directly incorporating a form 
of z dependence into the kinematic viscosity. But first, to better demonstrate that 
the zeroth order results are hardly affected by the choice of the z dependence in the 
kinematic viscosity, let us write down the result for the height of the (polytropic) 
standard alpha disk, where dv/dz = and 

v = ac s H. (2.35) 



If c s taken to be the zeroth order equatorial sound speed, c s o( r ) = c s o(r,0), one 
obtains the following zeroth order expressions for the kinematic and dynamic viscosity 
coefficients: 



Vo(r,z) = u p 



v/3 V^ 3/2 7 

h 2 (tf _ z 2fl* 



a 



5V15 



(2.36) 



(2.37) 



Notice that with this standard a-disk viscosity law, z/ depends solely on r and, 
therefore, r] inherits its z dependence entirely from the density, p (r,z), as given by 
eq. 

Now the disk semi-thickness, h(r), can be determined as a function of radius. To 



do this we observe that through lowest order in e, eq. ( 1.15 ) reads: 

dQ 



-h 



r] dz 



dr 



(2.38) 



where rh is the scaled constant from eq. Q2.34| ), j = r 1 / 2 is the (zeroth order) Keplerian 

i ii 

specific angular momentum, and j + = r + is the integration constant that arises from 



the no-torque boundary conditionP] that accompanies eq. (|1.15|) . Using eq. (|2.37 ) for 



3 If for some reason O was a function of both r and z we would need to use eq. ( |1.14| ). 
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r] and integrating over z, we are left with the following simple algebraic equation for 
the dimensionless disk height, h(r): 



h(r) 



Ar 



r 



1/6 



with 



An 



m / 80 5 
~a 1 37V 3 



1/6 



(2.39) 



This result implies that h(r)/r — > An =constant as r — > oo, and that the disk remains 
thin (i.e. H(r)/r ~ e) for all radii, provided, of course, that a is not too small. In 
addition, we see that as r — > r + , h — > 0. This is a consequence of our using j and 
dfl /dr in eq. ( |2.38| ), despite the requirement that dfl/dr — > as j (not j ) ^ J+- As 
we will see in Section 4, this is in fact a signal that our use of eq. ( p. 38; ) to determine 
h is not appropriate in the neighborhood of r + . 

It should be pointed out that the results we have so far obtained in this and 
the previous subsection are well known (c/. Shakura & Sunyaev 1973, Hoshi 1977, 
Paczyhski 1991). Novel developments only become apparent when eqs. (|1.21| )- (|1.24j ) 
are solved to higher order. But, as already remarked, eqs. ( |2.36|) and ([2.371) cannot 
be used to consistently extend the results obtained so far to higher order, it is first 
necessary to slightly modify the viscosity prescription. 



2.4. Lowest order results with height-dependent kinematic viscosity. 



In the original paper by Shakura Sz Sunyaev (1973) the dominant component of 
the viscous stress tensor is presumed to have the following form: 

o t4> ~ r]r(dVL/dr) a P. (2.40) 

Since we know that in the outer disk r{dVL/dr) ~ — f^, we use eq. ( |2.40| ) as a guide to 
make the following assumption regarding the z dependence of the viscosity: 

t ^ a 



c s (r, z) 
(1 + 1/n) 



(2.41) 
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consistent with r](r,z) ~ aP(r, z)/Qk where P ~ pc 2 for a polytrope. Our new 
expression for z/(r, z) reduces to the original formulation of eq. (|2.35|) in the midplane 
of the disk (z = 0). We justify our choice of v by noting that if the turbulent speed, 
Vturb, is bounded from above by the local sound speed, we must expect v tur b to vary 
considerably with height, since for a polytrope c s — > as z — > ±h. 



By using eqs. fl2.30|) -( f2.32|) & ( p. 41] ), we derive new forms for the zeroth order 
kinematic and dynamic viscosity coefficients: 



2a 

Un(r, z) = — 



Vo(r, z) = u (r, z)p (r,z) 



h 2 -z 2 

2a 
75^/5 



{h 2 



.2^ 5 /2' 



(2.42) 



;2.43) 



Comparison of eq. ( 2.43|) with eq. ( J2.37] ) shows that 7] is now one power higher in 
(ti 2 — z 2 ) and, as shown by Kita 1995, it is this difference that is the key to suppressing 
the divergence of v r and v z at the disk surface. 



We must first examine the impact that eqs. ( |2.42| )-( f2.43| ) may have on all other 
previously derived quantities. Since none of the results in § |2]^ depend in any way on 
v or 7], we know they are unaffected. The only change is in the value of the height of 
the disk (but not its functional form), which is increased by a factor 3 1//4 : 



h(r) 



A(l-W^ 
r 



1/6 



with A 



m /16(5) 
a \ 



3/2 \ ] V6 



7T 



, m 
1.96 ( — 

a 



1/6 



;2.44) 



2.5. Second order disk equations. 



To explore the differential rotation with respect to z and the nature of the 
velocity vector field in the accretion disk, we now consider only terms of 0(e 2 ) in 
eqs. ( |1.21| )-( |1.22| ). Bearing in mind that f2o = r _3//2 , Qi = 0, and uq = v o = v\ = 0, we 
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discover the following equations for Q 2 , U\, and v 2 : 

^ 3 z 2 3dc 2 1 9 / duA , n . rS 

-2^ 2 r = ----^ + --^ — j, (2.45) 

- — (rp u 1 ) + —( Po v 2 )=0. (2.47) 

Here, c s o, po and r] are given by eqs. fl2.30|) , (|2.42|) and ( |2.43|) , and h is known up to 
the integration constant r+. Unfortunately, the two viscous terms at the end of eqs. 
Q2.45|) & ( [2.46D complicate things by coupling the two equations together. To simplify 
the equations, most authors assume a priori that Q and v r in the outer disk are both 
functions of only r. As we will show, this is justified only in the limit a « e. We 
prefer to keep all terms so that we can obtain a solution valid through order e 3 which 
is consistent for all values of a. 



2.6. Complete analytical solution for 2 , u±, and v 2 . 

To solve for Q 2 and u\, we will make the ansatz that 
Ux(r,z) = fi(r)(h 2 — z 2 ) + /2(f), where /i(r) and f 2 (r) are as yet undetermined 
functions of r. A heuristic justification for this choice is that if the equations are 
decoupled by neglecting the z derivatives, as in the olden approach common in the 

literature, the solution for the lowest order corrections to Keplerian motion are 

2 



tt 2 \ old ( 



3 h 

4 \ r 



:| _ 2 dlnh 
3 V d In r 



'2.48) 



and Ui\ old (r, z) = g\{r){h 2 — z 2 ) + g 2 (r), where 



11 

in (y 



.5/2 



and 



92{r) 



-2a 



hr 

r 5/2 



din 

cHnr 



(2.49) 
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Given the nature of the coupling term, rjo(dui/dz), in eq. ( 2.45 ), it is now possible to 
formulate f2 2 in terms of /i(r) and /2(f), i-e. 



n 2 (r, z) = n 2 \° ld (r) + (^-) (h 2 - 6z 2 ) 



(2.50) 



We then solve for fi(r) and f2(r) by substituting eq. Q2.50| ) directly into the term 
involving r] (dQ2/ 'dz) that appears in eq. fl2.46p . This results in the following forms for 
the radial functions: 



fi(r) 



9i{r) 



and 



Mr) 



32 
— < 
15 



-ft' 



9i(r)h 2 \ 

1 I 64 ~2 +92{r), 



:a 2 



(2-51) 



where the functions g\{r) and g2{r) are defined by eq. ( |2.49| ). In this way we finally 
obtain complete solutions for Q2 and u\ in closed analytical form: 



n 



2 r 



uur, z) = —ft 



3 1 fdlnh\ 2 2l I, „z 2 ' 



dlnr 

where A is a constant that depends on the parameter a and is given by 

64 



(2.52) 
(2.53) 



5 / 



1 H ft' 

25 



(2.54) 



Note that Q,2(r,z) exhibits differential rotation with respect to the vertical 
coordinate; a feature which was also observed in the numerical solution by Kley & 
Lin (1992). This is because by including the viscous term of r]o(dui/dz) in the radial 
momentum equation we are no longer solving for Q under an assumption of strict 
radial hydrostatic equilibrium. 

We also observe that in the limit of ft 1, so that a 2 is a vanishingly small 
quantity, eqs. ( |2.52| )-( |2~5*3]) reduce identically to the solutions for Vt 2 and U\ of the 
equations of motion without the coupling terms. This suggests that neglecting the 
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viscous coupling terms is justified, so long as a is not too large, i.e. a < e. However, 
when a ~ 10 _1 to 1, the effects of including the 0(e 2 ) viscous terms becomes readily 
apparent. 

Note that Ui(r, z), being quadratic in z, is an even function with respect to 
reflections through the midplane (z = 0). It is also clear, upon reinstating the 
appropriate scale factors, that v r /v^k ~ H 2 /r 2 ~ e 2 in the outer disk, where our use 



of eq. (|2.44j) for the disk surface is known to be valid. Since (d\nh/d\nr) diverges at 
r + , we observe that lim r ^ r+ Ui(r, z) — > — oo, i.e. we have not cured the well known 
divergence of the Shakura-Sunyaev disk at the zero torque radius: as r — > r+, h — > 0, 
p — > and, to preserve M =const, v r — > oo. However, for all radii r > r+, u is finite 
everywhere and on the surface of the disk has the finite nonzero value given by 



u\(r, ±h) = —a 



r 5/2 



' dlnh 



352« 2 



d\nr J 3(25 + 64a 2 ) 



< 0. 



;2.55) 



With the knowledge that U\ remains finite at the surface of the disk, we can now 
pursue, with confidence, a solution for the lowest non-vanishing order of v z , i.e. by 
means of eq. (|2.47|) , the equation of continuity up through second order in e. If we 
use eq. Q2.31J) for p an d eq. (|2.53|) for %, then after the requisite differentiation (with 
respect to r) and integration (with respect to z), we can determine V2(r, z) up to an 
unknown function of the disk radius. Yet since v z is an odd function of z and therefore 
v 2 = in the equatorial plane for all r, this implies that the unknown function of r 
must be zero everywhere in the disk. In this manner, we find the following unique 
solution for V2- 



v 2 (r : z) = -a - 



.5/2 



" A I 1 -!?. 



32a 2 A fdlnh^ 



15 



dlnr 



' dlnh" 
dlnr 



(2.56) 



where A is the same constant that appears in eq. ( |2.54| ). 



We immediately notice several important features of the solution. First, as with 
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fl 2 (r) and ui(r, z), lim r ^ r+ t> 2 (r, z) — > — oo because of its dependence on (dlnh/dlnr) 2 . 
Secondly, we see that \v z /v r \ ~ H/r ~ e, as was expected for a standard thin disk in 
vertical hydrostatic equilibrium. Finally, v 2 , like U\, is finite along the surface of the 
disk for all r, i.e. 



v 2 (r, ±h) = =F« 



h 3 \ dlnh 



" 7 / 2 I din 1 , 



'dlnh\ 352a 2 



dlnr J 3(25 + 64a 2 



(2.57) 



We now check our solutions to see if they concur with what we expect for a standard 
thin disk. First, it is possible to show that eq. fl2.53| ) satisfies eq. ( |2.34| ) for the 



vertically-integrated mass flux, rh, even with the unexpected dependence on a 2 . 
Second, a quick glance reveals that our new solutions for Q 2 , v>i, and v 2 have the 
necessary parities with respect to reflection through the (z = 0) midplane, i.e. even, 
even, and odd, respectively. Third, if we replace the appropriate dimensional units for 
each quantity, we see that (v<j> — v^/v^j, and v r /v^k are both of C(e 2 ), while v z jv^ 
is 0(e 3 ); all of which is entirely in complete agreement with our assumption that the 
flow in the outer portions of the disk is predominantly Keplerian. 

Finally, since h(r) is the semi-thickness of the disk, we must have v z /v r = ±(dh/dr) 
for all points (r, ±ti) on the surface of the disk, as otherwise the geometrical disk 
surface cannot be in a steady-state. Comparison of eqs. (|2.55|) & ( [2. 571) shows that 



this constraint is indeed satisfied for the lowest non- vanishing orders of v r and v z : 

± v 2 (r,±h) = (h\( dlnh \ = dh 
U\(r, ±/i) \r J \d In r J dr 

We stress that the new solutions for Q 2 , u\, and v 2 fully comply with the above list 
of conditions for all possible values of a, provided that the disk remains geometrically 
thin. The remaining properties of eqs. ( |2.52| )-( |2~5"7|) for Q 2 , u\, and v 2 will be discussed 
in detail in the following section. 
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3. Detailed discussion of analytical results. 



In this section we examine the detailed properties of Vt 2 (eq. ||2.52 |), and in §f| 



those of the remaining components of the vector field, U\, and v 2 , eqs. ( ]2.53| )-( |2"37j ). 
Since by assumption ii) (§§ [!. 1| , |2.1| ) the next order corrections vanish everywhere in 
the disk, ^3 = ^2 = ^3 = 0, our equations are valid up to corrections of 0(z 4 /r 4 ) for 
Q and v r , and 0(z 5 /r 5 ) for v z . In other words, the results we describe are valid in 
the outer disk (r > 1.4r + ) up to relative corrections of e 2 , e.g. for H/R « 0.1 the 
expressions are valid to 1%. 



In § |3.1| , we study the nature of the z dependence of £l 2 (r, z). In § |3]3| we also 
discuss the interpretation of the apparent singularity in Q 2 at r = r+. In § |], we look 
at the velocity vector field in the outer disk, paying special attention to the sign of v r 
and v z . We find that beyond a certain radius there is significant mass outflow near 
the equatorial plane for a wide range of a's. We carefully analyze all of its important 
characteristics, including, in § [Q| , the fraction of the total mass flow carried out to 
infinity by this phenomenon. 



3.1. Analysis of the lowest-order correction to fi. 

We found that the angular velocity of material in the disk, up to and including 
terms of 0(e 2 ) is given by 

tt = tt jl + e 2 

Eq. (|3.59|) for Q is valid up to corrections of (9(z 4 /r 4 ) for a geometrically thin disk. 

As pointed out in Ch. |2] Q is clearly an even function of z and is Keplerian 
up through first order in e. In addition, since lim r .^ 00 ((iln/i/(ilnr) = 1 and 




h 2 



(3.59) 
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lim r _ >00 (/i/r) = A from eq. ( 2.44)) , we find that 



lim ( 77- - 1 



2 \ 2 



e z A 



4 15 1 h 2 , 



(3.60) 



Note that this limit is negative for all z because < 2a 2 A/15 < 0.082 over the range 
< a < 1. Also observe that in the case of a <C 1, eq. ( |3.60|) reduces to a negative 
constant — e 2 A 2 /4. Thus, in the limit of inviscid flow (a — * 0), f2 is constant on 
cylinders (but subkeplerian for r >> r + ). Finally, because h 2 (d In h/d In r) — > 00 as 
r — > r + , we know that Q, as given by eq. ( |3.59| ), diverges at r = r + . It is on these last 
two properties that we now focus our attention. 



3.2. Differential rotation of the angular velocity with respect to z. 

In this subsection, we discuss the z dependence of fl We define the quantity 
AQ = [Q(r, ±H) — Q(r, 0)] to be the difference between the angular velocity at the 
surface of the disk and the angular velocity in the equatorial plane. From eq. fl3.59| ). 
we derive the following expression for the fractional difference: 

(Afi/fio) = [0 2 (r, ±h) - n 2 (r, 0)] = -e 2 (j) Q« 2 a) . (3.61) 

Clearly, this fraction is negative for all radii. In Fig. 1, we plot log|Afi/f2o| m 
the disk with e = 0.01 and m — 1, for several values of the viscosity parameter: 
a = 0.01, 0.1, and 1.0. We observe that in the limit of r ^> r+, since lim r ^ 00 (/i/r) = A 
(see eq. ( |2.39| )), the fractional difference in Q tends to a negative constant given by 
lim r ^ oo (Afi/fi ) -> -e 2 (4/5)a 2 AA 2 (dashed lines). 

The amount by which the surfaces of constant Q deviate from upright cylinders 
is proportional to e 2 and therefore very small for a thin disk. The nearly constant 
behavior for (AQ/Qq) for r ^> r + is also interesting in its own right since it implies 
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that the total departure of the surfaces of constant f2 from the vertical does not 
significantly decrease (or increase) as one goes further out in the disk. Examination of 
eq. ( |3.61| ), also shows us that formally differential rotation vanishes as r — ► r+, since 
h — > at r = r + ; however, as Q = Qq + e 2 ^2 diverges to +00 at r = r + , this result of 
AQ — > really has no physical meaning. 



3.3. The singularity in fl(r, z) in the inner disk. 

In Fig. 2(a), we plot the second order correction to Q (i.e. e 2 ^) as a function of r 
in the (z = 0) midplane for three different values of a. In Fig. 2(b), we plot H = eh(r) 
for the same values a used in Fig. 2(a). For each choice of a, observe the singularity 
in Q 2 at r + . Clearly, as the presumed 0(e 2 ) correction provided by VL 2 grows without 
bound near r = r+, the assumptions underlying our perturbative expansion are invalid 
there. 

The reason for the divergence, at r = r + , of Q% , v r and v z can be traced back 
directly to the lowest order form of the vertically-integrated angular momentum 
equation ( f2.3£| ): 

™>Uo ~ 3+) = -r 3 



+h 

r] dz 

-h 



(3.62) 



dr 



The left side of eq. ( |3.62j ) is proportional to (j — j+) = jo (l — Jr+/rj and as 



such must vanish at r = r + . Since we are still assuming that Q w in the 
neighborhood of r + , if the right side is to also vanish there, then it is necessary for 
the height-integrated viscosity, J^rjodz, to go to zero at r = r+. However, in the 
a-disk prescription j^rjodz is proportional to (h/r) 6 , so this leads to the conclusion 
that h — > as r — > r + (Shakura & Sunyaev 1973). This functional form of h(r) causes 
a singularity in its derivative at r + , i.e. \im r ^ r+ (dh/dr) — > +00, and hence to a "cusp" 
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in the disk surface. 

The fault for all of this lies in the assumption that one can continue the 
perturbative expansion into the transition region, near the zero-torque point, r max , 
where Q reaches a maximum. While it is true that eq. ( |1.15| ) holds for all radii, 



including r = r + , the same cannot be said for its lowest order expansion in e, eq. ( |3.62| ), 
which was used in the derivation of h(r) and hence all other physical quantities. Since 
Q is maximal at r = r max , dVt/dr — > there, and eq. (|1.15|) is easily satisfied without 



requiring that t]odz vanish anywhere in the disk. However, eq. (|3.62|) , presumes 
that Q and j were approximately equal to their corresponding Keplerian values of fi 
and jo- This then forced the vertically-integrated viscosity (and hence h) to be zero at 
the radius, r + ; ultimately leading to the calculated singularities in fl, v r and v z . 

In the final analysis, we arrive at the conclusion that our perturbative expansion 
in e is not valid in the inner part of the disk, near r = r+, if we use the simple analytical 
expression for h(r) which appears in eq. ( |2.39| ). However, as they are written in terms 
of the disk surface, h(r), our solutions for Q 2 , u i, and v 2 in eqs. ( |2.52j )-( |2~56|) are more 



general than they first appear. Numerical work (Kita & Kluzniak 1997) shows that the 
radius of convergence for our expansions can be extended well into the boundary layer 
(r + ^ r) and all singularities completely disappear if only a more realistic (numerical) 
solution for h and (din h/dlnr) is used. In any event, our solutions for the vertical 
structure are certainly valid in the outer regions of the disk for r ^ 2r + , where we are 
fully justified in assuming Q w Q^, and hence in using eq. ( |3.62| ) to determine h. 
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4. The velocity vector field. 

In this section we will examine the behavior of the horizontal and vertical 
components of the fluid velocity in the accretion disk. In particular, we will concern 
ourselves with the sign of v r and v z , so that we can determine the direction of 
the accretion flow. Note that we adopt here the convention that accompanied 
eq. ( 1.12() , i- e - M > for accretion. Thus, if v r < the radial component of 



the flow is directed towards the central star. Likewise, v z < for z > signifies 
flow towards the z = midplane. Eqs. (p.53|) -( 2.57 ) for U\ and Vi directly give 



the unsealed velocity components as functions of r and z to lowest order in e, 
i.e. v r (r,z) = e 2 UiQk*R and v z (r, z) = e 3 t> 2 fifc*-R. 



4.1. The sign of v r in the outer disk: outflow vs. inflow. 



It is a result of our analysis, that at the surface of the accretion disk considered 
here, the fluid always flows in the general direction of the central object. Indeed, since 
2 > (32/15) Aa 2 for all possible a's in the range < a < 1, then v r < for all points 
on the disk surface because the term in eq. ( |2.53| ) proportional to (1 — z 2 /h 2 ) vanishes 
at z = ±/i. However, near the midplane beyond a certain radius, there may or may 
not be outflow depending on the chosen value of the parameter of viscosity, a. 



The radial velocity in the equatorial plane is: 



v r (r, 0) = — ae 2 £lk*R 



hT 



cHnr 



AH a' 

15 



(4.63) 



where (dlnh/dlnr) can be obtained from eq. (2.44) and Qk*R provides the physical 
units. For small radii, r ~ r+, evaluation of eq. ( [4.63| ) shows that the term involving 
(d In h/d In r) dominates the bracketed expression and thus v r (r, 0) < 0, i.e. we have 
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inflow. However, for large radii, r ^> r + , we see that since lim r ^, 00 (dlnh/dlnr) = 1, 
the sign of v r (r, 0) depends ultimately on whether or not 2 > A(l + 32a 2 /15). 

If, for example, a = 0, then A reduces to 11/5 and the bracketed portion of 
eq. ( |4.63|) is negative, so that v r (r, 0) > for r 3> r + ; indicating the existence of 



outflow in the equatorial plane far from the central accretoif]. On the other extreme, 
if a = 1 then A « 0.618 and A[l + 32a 2 / 15 ] ~ 1-94 < 2 and so v r (r, 0) < for r > r+. 
signaling that the equatorial flow is directed inwards towards the central star for all 



radii. A more precise analysis of eq. ( }4.63|) allows us to determine the critical value 



«cr = y 15/32 ~ 0.685 above which there is no backflow. This is a rather substantial 
value more than ten times as large as the critical value estimated by Kley & Lin 
(1992). 

We now analyze quantitatively where the direction of the flow changes sign when 
a < a cr . To this end we denote the stagnation radius, r stag , as the radius for which 
v r = in the equatorial plane. Eq. ( |4.63| ) with a < a cr allows us to calculate this 
stagnation radius as a function of a: 

2 

(4.64) 



stag 



a 



1 + 6 (A(l + ffa 2 ) - 2 



[6 (A(l + go") - 2) 

Observe that as a — > a cr = J (15/32), A — ► 1 by eq. ( |2.54j ), and r stag — > oo, as 



expected. This behavior is clearly visible in Fig 3, where we have plotted log (r sta9 /r + ) 
vs. a. In the opposite limit a — > 0, r stag /r + — > 121/36 ~ 3.36. For a < 10 -1 , we see 
that r stag w 3.5r + , this suggests that there is mass outflow in the midplane throughout 
most of the disk for realistic values of a. 



4 To our knowledge, equatorial outflow in accretion disks was first discovered by 
Urpin (1984). 
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Of course, for a < a cr , the region of outflow is not restricted to the midplane 
(z = 0) . Indeed, inspection of eq. ( [2.53| ) shows that v r (r, z) > for a range of z above 
and below the equator for r > r stag - We, therefore, define the "vertical flow surface," 
Zvert{r, a), as the surface on which v r = 0, implying that the flow is moving vertically 
there and hence the name. Clearly, z vert is well-defined for r > r stag and intersects 
the equatorial plane at r = r stag . Thus we can expect there to be outflow in the disk 
for all points contained in the domain of r > r stag {(y) and z < z ver t(r,a), whenever 
a < a cr . For a > a cr the surface of vertical flow disappears entirely. 



Solving for z vert via eq. (|2.53|) for m, we obtain 



,2 *)-U d ^l (4.65) 



/ z vert (r,a) \ _ 

V K r ) J ~ \ V ' 15^ J A \d\nr 
In Fig. 4, we show \og(Z vert /r + ) (solid curves) and \og(H/r + ) (dashed curves) for two 
different values^] of the parameter a: a = 0.1 & a = 0.6, where Z vert = ez vert and 
H = eh. Observe that in Fig. 4, since a = 0.6 is closer to a cr ~ 0.685, the outflow 
region, bounded from above by z vert (r,a), is beginning to be strongly suppressed and 
does not even start in the equatorial plane until r ^ 60r + . However, for a = 0.1 in 
Fig. 4, the region of backflow is quite extensive and occupies nearly 30% of the disk 
for r > 10r_i_. 



4.2. The sign of v z and its impact on mass outflow. 

We turn now our attention to the vertical velocity component, v z . A study 
of eq. ( |2.56|) reveals that ±v z < for all points on the disk surface at z = ±h(r), 
respectively, signifying that the flow there is directed towards the equatorial plane. 

5 H is also affected by the choice of a, since h oc (m/a) 1 / 6 . 
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In addition, eq. Q2.56 ) leads us to the conclusion that for a < a cr , besides being zero 
for all radii in the (z = 0) midplane because of its odd parity, v z also vanishes on a 
new and different "horizontal flow surface", Zhor(r,cx), on which the flow is directed 
horizontally. 

We determine the following functional form for this surface: 
/ z hor (r,a)\ 



32a 2 dlnh\ 2 ( dhiti" 2 



h{r) J \\ 15 d\nr J A\d\nr 



Comparison of this result with eq. (|4.65 ) reveals some definite similarities between 
the two flow surfaces, including the property that Zh or vanishes entirely in the limit 
a — > a cr . We also observe that z hor , the bounding surface across which v z changes 
sign, is contained entirely within the vertical flow surface, z vert , at which v r — 0. 
Furthermore, as r — > oo, Zh or asymptotically approaches z ver t from below for all values 
of a. Thus for a < a cr , there is a nested volume contained within the domain of 
backflow {y r > 0), within which the flow is directed away from the midplane,i.e. 
v z /z > 0. 

Figs. 5(a)-(b), illustrate this phenomenon for a = 0.1 and a = 0.6 respectively, 
where we have zoomed in to that fraction of the disk involved with outflow. Observe 
that at smaller radii, in both cases region A, the area within which v z /z > 0, is much 
smaller than the corresponding region B within which v r > 0; though their boundaries 
do converge to another at infinity, since lim r .^ 00 (<iln/i/(ilnr) — > 1. 



Eqs. ( p.53| )- (|2.56|) , in conjunction with eqs. Q4.65| )- (|4.66|) , also allow us to draw 



some rather interesting conclusions regarding the flow geometry in the outer disk. 
First, any material that enters region B must cross over into region A at some time 
in the future. To see this, note that v z /z is still negative for all points between the 
vertical flow surface and the horizontal flow surface (l^orl < \z\ < \z ver t\) an d thus 
any fluid element in region B (where v r > 0), given enough time, must eventually 
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cross the inner horizontal flow surface, Zh or , and into region A. It is clear that some of 
the material in region C in Figs. 5(a)-(b) must not enter region B, instead a sizeable 
fraction of the flow that exists above z vert must continue on ahead of the stagnation 
radius, into region D and onto the central star. This is, of course, the accreted flow 
that comprises M. 

To better visualize our previous remark, we call the reader's attention to 
Figs. 6(a)-(b), where we plot the direction of the velocity vector field, i.e. the unit 
vectors formed from the components v r and v z , in the outer disk for two different 
values of viscosity, both with a < a cr . As a reference, we have also overlaid the 
vertical flow surface, z vert . In Fig. 6(c), we show the directional flow pattern for the 
case of a > a cr (here a — 1), and we now find that only inflow persists throughout the 
entire disk. 

Finally, note that there are no cells of meridional circulation in the disk. In the 
region B flow is always directed towards the surface which bounds region A, while in 
the region A flow is always directed away from the equator. The circulating pattern 
of flow found by Kley & Lin (1992) close to the edges of their computational domain 
must be an artifact of the boundary conditions they imposed. In fact, calculation of 
the stream lines shows that material which enters region A can never intersect the 
v z — horizontal flow surface again, because v r > and v z — > there. Coupled with 
our original finding that any fluid which enters region B must also pass eventually into 
region A, this implies that material once in the backflow region B can never return to 
the inflow portion C of the accretion disk. Therefore, if a < a cr , the backflow must 
continue all the way out to "infinity," where the disk terminates. All of this can be 
easily verified by looking at Figs. 6(a)-(b). 
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4.3. The mass fraction of outflow relative to inflow. 



The last question that we address relates to the backflow. What is the mass 
outflow rate in comparison to the M accreted by the star? It is clear that the answer 
to this question must depend on the parameter a, since we know that there is no 
backflow for a > a cr ~ 0.685. For this reason we now define T(r, a) to be the ratio of 
the mass rate flowing outwards through a cylinder of radius r > r stag to the net mass 
accretion rate, M. 

To evaluate the functional form of T(r,a), we simply integrate the radial mass 
flux, pv r , over the surface of a cylinder with height 2z ver t centered on the midplane. 



r(r, a) 



4nr f^(r,«) pQUidz 



1 






2r [ 


rh 


Jo 



rzvert(r,a) 



4vrr J /l{r) poUidz 

where fn is the scaled mass accretion rate defined in eq. ( |2.34| ). 



pQU\dz 



(4.67) 



(4.68) 



Using eq. ( |4.65| ) for z vert (r, a), eq. ( |2.31| ) for p , and eq. ( |2.53| ) for Ui(r, z) we find 
that the mass outflow fraction is 

T(r a) - ( 2aA ) ( k ) 6 
\5Vbm J \rj 

where 7* = 7*(r, a) = z vert /h and G\ and G 2 are rather complicated functions of 7* 
given by: 

" 2 \ 5/2 dz 5 . 1 , ^ 5 



+ 



(Vi-(r) 2 ) 3 + ^(Vi-(r) : 



(4.69) 



+^7* (\/i-(r) 2 ) 3 . 



(4.70) 
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The limit r 00(a) = Hindoo T(r, a) is of particular interest, as it represents the ratio of 
the net outflow to net inflow in the disk. 

Analysis of eqs. ( |4.68; )- (|4.70| ) provides us with the following form for T^: 

32A 



7T 



Gi(tSo) + ((T^) 2 - l) G 2 tfJ\ , (4.71) 



where the quantity 7^ is well defined only for a < a cr and given by 



/ 32 2 

TSo = J™ 7.(r, «) = Wl + -« 2 - x . (4.72) 



The resultant dependence on a for the mass outflow fraction as r — > 00 is plotted 
in Fig. 7. We see that r oo (a) tends smoothly to zero as a — > a cr . However, for a wide 
range of a, the fraction of the mass flow that is contained in the outflow region, near 
the midplane is quite significant: r 00(a) ~ 0.4 for < a ^ 0.1 and 0.1 < T^a) ^ 0.35 
for 0.1 < a < 0.5. Thus, for a < 0.1, the total mass rate, M to t = M out + M, being fed 
into the disk at the outer edge is ~ 1AM. 

All in all, we are led to the startling conclusion (Urpin 1984) that, for small values 
of a, there is backflow in the disk transporting fluid outwards to the outer boundary 
of the accretion disk. This should have serious repercussions for mass transfer at the 
outer edge of the disk. 
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6. Final remarks. 

By performing a systematic expansion (pioneered by Regev 1983) of the disk equations 
of motion we were able to find a closed solution for the velocity field and the disk 
structure, valid everywhere outside, say, 1.1 times the zero torque radius. We showed 
how, for all but very large values of viscosity, the accretion flow turns around and 
feeds a backflow (first discovered by Urpin in 1984) in the equatorial plane of the disk. 
This backflow has now also been seen in a number of numerical simulations and must 
be considered a general feature of accretion in a geometrically thin disk, and possibly 
also in quasi-spherical flows. 

We note, that if the flow discussed here were advective, some of the gravitational 
energy released by the flow close to the central gravitating body would have been 
carried by the flow to larger radii before it is radiated. An urgent topic of investigation 
should be whether solutions with backflows, similar to the one presented here, may be 
present in an advection dominated flow. Should such solutions exist, conclusions (e.g. 
Narayan 1996) that the apparent deficit of emission in the inner region of accretion 
disks of some X-ray "novae" necessarily implies the presence of a space-time horizon 
would have to be treated with caution. 
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Fig. 1. — The fractional difference between the angular velocity at the disk surface 
and the angular velocity in the equatorial plane relative to Keplerian as a function 
of the radius, for e = 0.01 and a =0.01, 0.1, 1.0. The dashed lines correspond to 
linv^oolog |Af2/f2 |. Note that in all the figures we arbitrarily choose m — 1. 



Fig. 2(a). — Plot of the second order correction to Vt vs. radius for e = 0.01 and 
a = 0.01, 0.1, 0.5. Note the change in sign for r « 1.4r + for all three values of a. 



Fig. 2(b). — The surface of the disk, H = eh(r), for the identical set of parameters 
used in Fig. 2(a). Observe that for all three values of a, the singularity in Q in Fig. 
2(a) at r = r + corresponds to H — > 0. 
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Fig. 3. — The stagnation radius vs. a. Observe that for a < 0.2, r stag < 4.0r + (tending 
to 121/36 « 3.36 for a = 0), but diverges rapidly asa-* a cr = ^(15/32) « 0.685. 



Fig. 4. — Logscale plot of Z vert = ez vert (solid lines) and H = eh (dashed lines) for 
a — 0.1 and a = 0.6, with e = 0.01 in both cases. Notice how much smaller the vertical 
flow surface is relative to the disk surface when a = 0.6. 



Fig. 5(a). — Plot of Z vert = ez vert (solid lines), Zh or = e-^W (dashed lines), and H = eh 
(dot-dashed lines) for a = 0.1 and e = 0.01. Region A represents the zone where v r > 
and v z /z > 0, B labels the area where v r > and v z /z < 0, and C marks the portion 
where v r < and v z < 0. The arrows labeled D indicate the region where the fluid 
flows only towards the central object. 
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Fig. 5(b). — Same as Fig. 5(a) except now a = 0.6. Note the change in scale and 
observe how much smaller the two flow surfaces are relative to the disk surface when 
a = 0.6. 



Fig. 6(a). — The unit velocity vector field, i.e. (v r /\v r \,iT z /\v z \) for e = 0.1 and a = 0.1. 
The dashed curve represents the disk surface and the solid curve the vertical flow surface 
where v r — 0. 



Fig. 6(b). — Same as Fig. 6(a) except now a = 0.6. Note that because the outflow 
region has been greatly reduced and pushed out to r > 62.5r + , we have zoomed in to 
the equatorial region near the stagnation radius. 



Fig. 6(c). — Same as Figs. 6(a)-(b) except now a = 1.0. There is no longer any outflow 
in the disk as a > a cr ~ 0.685. 



39 



Fig. 7. — The ratio of the mass flow rate of material moving outwards (through a 
cylinder with r = oo) to the net accreted mass flow rate, M, vs. a. Observe that 
r 00(a) vanishes as a — > a cr rs 0.685, but is about 40% for a < 0.1. 



